library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores=3)

# model without candidate random effects
model1 <- "
data{
int N;
int K;

int P;
int<lower=1, upper=P> Pty[N];

vector<lower=0, upper=100>[N] y;

matrix[N,K] X;
}

parameters{
real alpha;
vector[K] beta;
vector[P] nu_p;

real<lower=0, upper=100> sigma_y;
real<lower=0, upper=100> sigma_p;
}

model{
sigma_y ~ uniform(0, 100);
sigma_p ~ uniform(0, 100);

y ~ normal(alpha + X*beta + nu_p[Pty], sigma_y);

nu_p ~ normal(0, sigma_p);
}
"

stanModel1 <- stan_model(model_code=model1)


# model with candidate random effects
model2 <- "
data{
int N;
int K;

int C;
int<lower=1, upper=C> Cand[N];

int P;
int<lower=1, upper=P> Pty[N];

vector<lower=0, upper=100>[N] y;

matrix[N,K] X;
}

parameters{
real alpha;
vector[K] beta;
vector[C] nu_c;
vector[P] nu_p;

real<lower=0, upper=100> sigma_y;
real<lower=0, upper=100> sigma_c;
real<lower=0, upper=100> sigma_p;
}

model{
sigma_y ~ uniform(0, 100);
sigma_c ~ uniform(0, 100);
sigma_p ~ uniform(0, 100);

y ~ normal(alpha + X*beta + nu_c[Cand] + nu_p[Pty], sigma_y);

nu_c ~ normal(0, sigma_c);
nu_p ~ normal(0, sigma_p);
}
"

stanModel2 <- stan_model(model_code=model2)




load("JHRED2014_New.RData")




###
### sntv vote share (models 1 to 6 in Table 1)
###

# election-year matrix
year_list <- sort(unique(sntv$year))
year.mat.sntv <- matrix(0, nrow=nrow(sntv), ncol=length(year_list)-1)

for(i in 2:length(year_list)){
  year.mat.sntv[sntv$year==year_list[i],i-1] <- 1
}

# explanatory variable matrix
# kakusu = name complexity (# of strokes)
# kakusu.ave = average name complexity (# of strokes/name length)
x.mat.sntv <- cbind(sntv$kakusu, sntv$length, log(sntv$ku_ncand), sntv$ku_m, 
                    year.mat.sntv)
x.mat2.sntv <- cbind(sntv$kakusu.ave, log(sntv$ku_ncand), sntv$ku_m, 
                     year.mat.sntv)

sntv.dat1 <- list(N=nrow(sntv),
                  K=ncol(x.mat.sntv),
                  y=sntv$voteshare,
                  X=x.mat.sntv,
                  P=length(unique(sntv$ptyID)),
                  Pty=sntv$ptyID)

sntv.all1 <- sampling(stanModel1, sntv.dat1, iter=3000, warmup=1000, 
                      thin=1, chains=3, seed=1234)

					  
sntv.dat2 <- list(N=nrow(sntv),
                  K=ncol(x.mat.sntv),
                  y=sntv$voteshare,
                  X=x.mat.sntv,
                  C=length(unique(sntv$candID)),
                  Cand=sntv$candID,
                  P=length(unique(sntv$ptyID)),
                  Pty=sntv$ptyID)

sntv.all2 <- sampling(stanModel2, sntv.dat2, iter=3000, warmup=1000, 
                      thin=1, chains=3, seed=1234)

					  
sntv.dat3 <- list(N=nrow(sntv),
                  K=ncol(x.mat2.sntv),
                  y=sntv$voteshare,
                  X=x.mat2.sntv,
                  P=length(unique(sntv$ptyID)),
                  Pty=sntv$ptyID)

sntv.all3 <- sampling(stanModel1, sntv.dat3, iter=3000, warmup=1000, 
                      thin=1, chains=3, seed=1234)

					  
sntv.dat4 <- list(N=nrow(sntv),
                  K=ncol(x.mat2.sntv),
                  y=sntv$voteshare,
                  X=x.mat2.sntv,
                  C=length(unique(sntv$candID)),
                  Cand=sntv$candID,
                  P=length(unique(sntv$ptyID)),
                  Pty=sntv$ptyID)

sntv.all4 <- sampling(stanModel2, sntv.dat4, iter=3000, warmup=1000, 
                      thin=5, chains=2, seed=1234)

					  
## sntv first-time candidates only
sntv.first <- sntv[sntv$totcruns==1,]

# election year matrix
year_list <- sort(unique(sntv.first$year))
year.mat.sntvfirst <- matrix(0, nrow=nrow(sntv.first), ncol=length(year_list)-1)

for(i in 2:length(year_list)){
  year.mat.sntvfirst[sntv.first$year==year_list[i],i-1] <- 1
}

# explanatory variable matrix
x.mat.sntvfirst <- cbind(sntv.first$kakusu, sntv.first$length, log(sntv.first$ku_ncand), 
                         sntv.first$ku_m, year.mat.sntvfirst)
x.mat2.sntvfirst <- cbind(sntv.first$kakusu.ave, log(sntv.first$ku_ncand), sntv.first$ku_m,
                          year.mat.sntvfirst)

sntv.dat5 <- list(N=nrow(sntv.first),
                  K=ncol(x.mat.sntvfirst),
                  y=sntv.first$voteshare,
                  X=x.mat.sntvfirst,
                  P=length(unique(sntv.first$ptyIDfirst)),
                  Pty=sntv.first$ptyIDfirst)

sntv.first1 <- sampling(stanModel1, sntv.dat5, iter=3000, warmup=1000, 
                        thin=1, chains=3, seed=1234)

						
sntv.dat6 <- list(N=nrow(sntv.first),
                  K=ncol(x.mat2.sntvfirst),
                  y=sntv.first$voteshare,
                  X=x.mat2.sntvfirst,
                  P=length(unique(sntv.first$ptyIDfirst)),
                  Pty=sntv.first$ptyIDfirst)

sntv.first2 <- sampling(stanModel1, sntv.dat6, iter=3000, warmup=1000, 
                        thin=1, chains=3, seed=1234)

#save(sntv.all1, sntv.all2, sntv.all3, sntv.all4,
#     sntv.first1, sntv.first2, file="StanFits_SNTVComplexity.RData")




###
### smdp vote share (models 1 to 6 in Table 2)
###

# election-year matrix
year_list <- sort(unique(smdp$year))
year.mat.smdp <- matrix(0, nrow=nrow(smdp), ncol=length(year_list)-1)

for(i in 2:length(year_list)){
  year.mat.smdp[smdp$year==year_list[i],i-1] <- 1
}

# explanatory variable matrix
x.mat.smdp <- cbind(smdp$kakusu, smdp$length, log(smdp$ku_ncand), year.mat.smdp)
x.mat2.smdp <- cbind(smdp$kakusu.ave, log(smdp$ku_ncand), year.mat.smdp)

smdp.dat1 <- list(N=nrow(smdp),
                  K=ncol(x.mat.smdp),
                  y=smdp$voteshare,
                  X=x.mat.smdp,
                  P=length(unique(smdp$ptyID)),
                  Pty=smdp$ptyID)

smdp.all1 <- sampling(stanModel1, smdp.dat1, iter=3000, warmup=1000, 
                      thin=1, chains=3, seed=1234)

					  
smdp.dat2 <- list(N=nrow(smdp),
                  K=ncol(x.mat.smdp),
                  y=smdp$voteshare,
                  X=x.mat.smdp,
                  C=length(unique(smdp$candID)),
                  Cand=smdp$candID,
                  P=length(unique(smdp$ptyID)),
                  Pty=smdp$ptyID)

smdp.all2 <- sampling(stanModel2, smdp.dat2, iter=3000, warmup=1000, 
                      thin=1, chains=3, seed=1234)

					  
smdp.dat3 <- list(N=nrow(smdp),
                  K=ncol(x.mat2.smdp),
                  y=smdp$voteshare,
                  X=x.mat2.smdp,
                  P=length(unique(smdp$ptyID)),
                  Pty=smdp$ptyID)

smdp.all3 <- sampling(stanModel1, smdp.dat3, iter=3000, warmup=1000, 
                      thin=1, chains=3, seed=1234)

					  
smdp.dat4 <- list(N=nrow(smdp),
                  K=ncol(x.mat2.smdp),
                  y=smdp$voteshare,
                  X=x.mat2.smdp,
                  C=length(unique(smdp$candID)),
                  Cand=smdp$candID,
                  P=length(unique(smdp$ptyID)),
                  Pty=smdp$ptyID)

smdp.all4 <- sampling(stanModel2, smdp.dat4, iter=3000, warmup=1000, 
                      thin=1, chains=3, seed=1234)

					  
## smdp first-time candidates only
smdp.first <- smdp[smdp$totcruns==1,]

# election year matrix
year_list <- sort(unique(smdp.first$year))
year.mat.smdpfirst <- matrix(0, nrow=nrow(smdp.first), ncol=length(year_list)-1)

for(i in 2:length(year_list)){
  year.mat.smdpfirst[smdp.first$year==year_list[i],i-1] <- 1
}

# explanatory variable matrix
x.mat.smdpfirst <- cbind(smdp.first$kakusu, smdp.first$length, log(smdp.first$ku_ncand), 
                         year.mat.smdpfirst)
x.mat2.smdpfirst <- cbind(smdp.first$kakusu.ave, log(smdp.first$ku_ncand),
                          year.mat.smdpfirst)

smdp.dat5 <- list(N=nrow(smdp.first),
                  K=ncol(x.mat.smdpfirst),
                  y=smdp.first$voteshare,
                  X=x.mat.smdpfirst,
                  P=length(unique(smdp.first$ptyIDfirst)),
                  Pty=smdp.first$ptyIDfirst)

smdp.first1 <- sampling(stanModel1, smdp.dat5, iter=3000, warmup=1000, 
                        thin=1, chains=3, seed=1234)

						
smdp.dat6 <- list(N=nrow(smdp.first),
                  K=ncol(x.mat2.smdpfirst),
                  y=smdp.first$voteshare,
                  X=x.mat2.smdpfirst,
                  P=length(unique(smdp.first$ptyIDfirst)),
                  Pty=smdp.first$ptyIDfirst)

smdp.first2 <- sampling(stanModel1, smdp.dat6, iter=3000, warmup=1000, 
                        thin=1, chains=3, seed=1234)

#save(smdp.all1, smdp.all2, smdp.all3, smdp.all4,
#     smdp.first1, smdp.first2, file="StanFits_SMDPComplexity.RData")




###
### sntv vote share with difficulty scores (models 1 to 6 in Table E.2)
###

# election-year matrix
year_list <- sort(unique(sntv$year))
year.mat.sntv <- matrix(0, nrow=nrow(sntv), ncol=length(year_list)-1)

for(i in 2:length(year_list)){
  year.mat.sntv[sntv$year==year_list[i],i-1] <- 1
}

# diff = difficulty score
# diff = averaege difficulty score (difficulty score/name length)
x.mat3.sntv <- cbind(sntv$diff, sntv$length, log(sntv$ku_ncand), sntv$ku_m, year.mat.sntv)
x.mat4.sntv <- cbind(sntv$diff.ave, log(sntv$ku_ncand), sntv$ku_m, year.mat.sntv)

sntv.dat7 <- list(N=nrow(sntv),
                  K=ncol(x.mat3.sntv),
                  y=sntv$voteshare,
                  X=x.mat3.sntv,
                  P=length(unique(sntv$ptyID)),
                  Pty=sntv$ptyID)

sntv.diff.all1 <- sampling(stanModel1, sntv.dat7, iter=3000, warmup=1000, 
                           thin=1, chains=3, seed=1234)

						   
sntv.dat8 <- list(N=nrow(sntv),
                  K=ncol(x.mat3.sntv),
                  y=sntv$voteshare,
                  X=x.mat3.sntv,
                  C=length(unique(sntv$candID)),
                  Cand=sntv$candID,
                  P=length(unique(sntv$ptyID)),
                  Pty=sntv$ptyID)

sntv.diff.all2 <- sampling(stanModel2, sntv.dat8, iter=3000, warmup=1000, 
                           thin=1, chains=3, seed=1234)

						   
sntv.dat9 <- list(N=nrow(sntv),
                  K=ncol(x.mat4.sntv),
                  y=sntv$voteshare,
                  X=x.mat4.sntv,
                  P=length(unique(sntv$ptyID)),
                  Pty=sntv$ptyID)

sntv.diff.all3 <- sampling(stanModel1, sntv.dat9, iter=3000, warmup=1000, 
                           thin=1, chains=3, seed=1234)

						   
sntv.dat10 <- list(N=nrow(sntv),
                   K=ncol(x.mat4.sntv),
                   y=sntv$voteshare,
                   X=x.mat4.sntv,
                   C=length(unique(sntv$candID)),
                   Cand=sntv$candID,
                   P=length(unique(sntv$ptyID)),
                   Pty=sntv$ptyID)

sntv.diff.all4 <- sampling(stanModel2, sntv.dat10, iter=3000, warmup=1000, 
                           thin=1, chains=3, seed=1234)

						   
# sntv first-time candidates only
sntv.first <- sntv[sntv$totcruns==1,]

# election year matrix
year_list <- sort(unique(sntv.first$year))
year.mat.sntvfirst <- matrix(0, nrow=nrow(sntv.first), ncol=length(year_list)-1)

for(i in 2:length(year_list)){
  year.mat.sntvfirst[sntv.first$year==year_list[i],i-1] <- 1
}

x.mat3.sntvfirst <- cbind(sntv.first$diff, sntv.first$length, log(sntv.first$ku_ncand), 
                          sntv.first$ku_m, year.mat.sntvfirst)
x.mat4.sntvfirst <- cbind(sntv.first$diff.ave, log(sntv.first$ku_ncand), sntv.first$ku_m, 
                          year.mat.sntvfirst)

sntv.dat11 <- list(N=nrow(sntv.first),
                   K=ncol(x.mat3.sntvfirst),
                   y=sntv.first$voteshare,
                   X=x.mat3.sntvfirst,
                   P=length(unique(sntv.first$ptyIDfirst)),
                   Pty=sntv.first$ptyIDfirst)

sntv.diff.first1 <- sampling(stanModel1, sntv.dat11, iter=3000, warmup=1000, 
                             thin=1, chains=3, seed=1234)

							 
sntv.dat12 <- list(N=nrow(sntv.first),
                   K=ncol(x.mat4.sntvfirst),
                   y=sntv.first$voteshare,
                   X=x.mat4.sntvfirst,
                   P=length(unique(sntv.first$ptyIDfirst)),
                   Pty=sntv.first$ptyIDfirst)

sntv.diff.first2 <- sampling(stanModel1, sntv.dat12, iter=3000, warmup=1000, 
                             thin=1, chains=3, seed=1234)

#save(sntv.diff.all1, sntv.diff.all2, sntv.diff.all3, sntv.diff.all4,
#     sntv.diff.first1, sntv.diff.first2, file="StanFits_SNTVDifficulty.RData")





###
### summary
###
library(rstan)
#load("StanFits_SNTVComplexity.RData")
#load("StanFits_SMDPComplexity.RData")
#load("StanFits_SNTVDifficulty.RData")

# Table 1
summary(sntv.all1, pars=c("beta[1]", "beta[2]", "beta[3]", "beta[4]", "sigma_p"))$summary[,c(1,3)]
summary(sntv.all2, pars=c("beta[1]", "beta[2]", "beta[3]", "beta[4]", "sigma_c", "sigma_p"))$summary[,c(1,3)]
summary(sntv.all3, pars=c("beta[1]", "beta[2]", "beta[3]", "sigma_p"))$summary[,c(1,3)]
summary(sntv.all4, pars=c("beta[1]", "beta[2]", "beta[3]", "sigma_c", "sigma_p"))$summary[,c(1,3)]
summary(sntv.first1, pars=c("beta[1]", "beta[2]", "beta[3]", "beta[4]", "sigma_p"))$summary[,c(1,3)]
summary(sntv.first2, pars=c("beta[1]", "beta[2]", "beta[3]", "sigma_p"))$summary[,c(1,3)]


# Table 2
summary(smdp.all1, pars=c("beta[1]", "beta[2]", "beta[3]", "sigma_p"))$summary[,c(1,3)]
summary(smdp.all2, pars=c("beta[1]", "beta[2]", "beta[3]", "sigma_c", "sigma_p"))$summary[,c(1,3)]
summary(smdp.all3, pars=c("beta[1]", "beta[2]", "sigma_p"))$summary[,c(1,3)]
summary(smdp.all4, pars=c("beta[1]", "beta[2]", "sigma_c", "sigma_p"))$summary[,c(1,3)]
summary(smdp.first1, pars=c("beta[1]", "beta[2]", "beta[3]", "sigma_p"))$summary[,c(1,3)]
summary(smdp.first2, pars=c("beta[1]", "beta[2]", "sigma_p"))$summary[,c(1,3)]


#Table E.2 in the appendix
summary(sntv.diff.all1, pars=c("beta[1]", "beta[2]", "beta[3]", "beta[4]", "sigma_p"))$summary[,c(1,3)]
summary(sntv.diff.all2, pars=c("beta[1]", "beta[2]", "beta[3]", "beta[4]", "sigma_c", "sigma_p"))$summary[,c(1,3)]
summary(sntv.diff.all3, pars=c("beta[1]", "beta[2]", "beta[3]", "sigma_p"))$summary[,c(1,3)]
summary(sntv.diff.all4, pars=c("beta[1]", "beta[2]", "beta[3]", "sigma_c", "sigma_p"))$summary[,c(1,3)]
summary(sntv.diff.first1, pars=c("beta[1]", "beta[2]", "beta[3]", "beta[4]", "sigma_p"))$summary[,c(1,3)]
summary(sntv.diff.first2, pars=c("beta[1]", "beta[2]", "beta[3]", "sigma_p"))$summary[,c(1,3)]